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ABSTRACT 

Fundamental aspects of spectral methods are introduced. Recent 
developments in spectral methods are reviewed with an emphasis on 
collocation techniques. Their applications to both compressible and 
incompressible flows, to viscous as well as inviscid flows, and also to 
chemically reacting flows are surveyed. The key role that these methods 
play in the simulation of stability, transition, and turbulence is 
brought out. A perspective is provided on some of the obstacles that 
prohibit a wider use of these methods, and how these obstacles are being 
overcome. 
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INTRODUCTION 


In certain areas of computational fluid dynamics spectral methods 
have become the prevailing numerical tool for large-scale calculations. 
This is certainly the case for such three-dimensional applications as 
direct simulation of homogeneous turbulence, computation of transition 
in shear flows, and global weather modeling. For many other applica- 
tions, such as heat transfer, boundary layers, reacting flows, compres- 
sible flows, and magnetohydrodynamics, spectral methods have proven to 
be a viable alternative to the traditional finite difference and finite 
element techniques. 

Spectral methods are characterized by the expansion of the solution 
in terms of global and, usually, orthogonal polynomials. Since the mid- 
nineteenth century this has been a standard analytical tool for linear, 
separable differential equations. Nonlinearities present considerable 
algebraic difficulties, even on a modern computer. These were 
surmounted effectively in the early 1970's, and only then did spectral 
methods become competitive with alternative algorithms. By the present 
time, however, spectral methods have been refined and extended to the 
point where many problems in fluid mechanics are only tractable by this 
technique. 

Numerical spectral methods for partial differential equations were 
originally developed by meteorologists. Though this approach was 
proposed by Blinova in 1943 and Haurwitz and Craig in 1952, the first 
numerical computations were conducted by Silberman (1954). The expense 
of computing nonlinear terms remained a severe drawback until Orszag 
(1969) and Eliasen, et al (1970) developed the transform methods that 
still form the backbone of many large-scale spectral computations. 
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These methods and others used in fluid mechanics prior to 1970 are 
now termed spectral Galerkin methods: the fundamental unknowns are the 
expansion coefficients and the equations for these are derived by the 
techniques used in classical analysis. The advent of computers made 
feasible an alternative discretization, termed the spectral collocation 
technique, in which the fundamental unknowns are the solution values at 
selected, collocation points and the series expansion is used solely for 
the purpose of approximating derivatives. This approach was proposed by 
Kreiss and Oliger (1971) and Orszag (1972). 

Many useful versions of spectral methods have been developed since 
1971 and especially during the 1980's. This review will discuss many of 
the recent innovations and will focus on the collocation technique since 
it is the version most readily applicable to nonlinear problems. We 
will survey applications to both compressible and incompressible flows, 
to viscous as well as inviscid flows, and also to chemically reacting 
flows. In the interests of brevity we shall not cover the applications 
to meteorology, magnetohydrodynamics, astrophysics, and other related 
fields. Moreover, we will restrict ourselves to the three-dimensional 
applications of well-established algorithms while discussing some two- 
and even one-dimensional applications of more novel spectral methods. 

Let us mention here some other articles for those interested in 
additional historical references, applications in other fields, and 
theoretical developments on the numerical analysis of spectral methods. 
The monograph by Gottlieb and Orszag (1977) describes the theory and 
applications developed prior to 1977. It will be referenced hereafter 
as GO. The following five years are covered in the proceedings edited 
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by Voigt, et al (1984), Fluid dynamical applications, especially 
multigrid techniques, are discussed by Zang and Hussaini (1985c). 
Compressible flow applications are covered by Hussaini, Salas, and Zang 
(1985). The role of spectral methods in meteorology is explained by 
Jarraud and Baede (1985). The book by Canuto, et al (1987) contains a 
detailed description of many spectral algorithms and presents an 
exhaustive discussion of the theoretical aspects of these numerical 
methods. It will be referenced hereafter as CHQZ. 


FUNDAMENTALS 

The motivation for the use of spectral methods in numerical 
calculations stems from the attractive approximation properties of 
orthogonal polynomial expansions. Suppose, for example, that a 
function u(x) is expanded in a truncated Chebyshev series on [-1,1]: 

N 

V*> “ £ V*> (1) 

n=0 

where T n (x) = cos(n arc cos x). The classical form of the expansion 
coefficients (or spectra) is 

1 

a = (2/c ) f u(x) T (x)(l - x 2 ) ^ 2 dx (2) 

n n J n v ' 

-1 

where Cq = 2, and c n = 1 for n > 1. The substitution x = cos 0 
converts this into a Fourier cosine series. A simple integration-by- 
parts argument (GO, Ch. 3) reveals that 
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h* 5 a n 0» as n «, for all p > 0 (3) 

provided that u is infinitely differentiable. Consequently, the 
approximation error decreases faster than algebraically. This rapid 
convergence is referred to as infinite order accuracy, exponential 
convergence, or spectral accuracy. Our primary concern in this review 
is on numerical methods for partial differential equations that exhibit 
spectral accuracy for infinitely differentiable solutions. 

The approximation just described is typical of spectral Galerkin 
methods. An alternative approximation, termed spectral collocation, is 
one of interpolation. It retains the expansion (1), but replaces the 
condition (2) for the expansion coefficients, with the condition 


u N (xj) = u(xj) (4) 

where xj are special, so-called collocation, points in [-1,1], For 
most problems, the optimal choice of these collocation points is 

Xj = cos(irj/N). ( 5 ) 

This choice of collocation points yields an extremely accurate 
approximation (CHQZ, Ch.,2) to the integral appearing in Eq. (2): 

N 

a „ - (2 /Nc„) 'Ey 1 u( yVV’ 

j=0 


( 6 ) 
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where, = c„ = 2 and c =1, otherwise. Whether (2) or (6) is used 
ON n 

for the expansion coefficients, the expansion (1) is differentiated 
analytically to form the approximations to whatever derivatives are 
required for the problem at hand. 

A graphical distinction between traditional approximations and 
spectral ones is provided in Figure 1 for the simple task of estimating 
the derivative of the function 1 + sin (2ttx + tt/4) on [-1,1] from the 
values of the function at a finite number of grid points. A finite 
difference or finite element method uses local information to estimate 
derivatives whereas a spectral method uses global information. In this 
figure a second-order (central) finite difference method is compared 
with a Chebyshev spectral collocation method. The finite difference 
approximation estimates the derivative at, say, x = 0, from the parabola 
which interpolates the function at x = 0 and the two adjacent grid 
points. A separate parabola is used at each grid point. The spectral 
approximation, on the other hand, uses all the available information 
about the function. If there are N + 1 grid points, then the 
interpolating polynomial from which the derivative is extracted has 
degree N and the same polynomial is used for all the grid points. 
Note that the local method produces a second-order accurate derivative, 
with the error decreasing as l/N^, whereas the error from the global 
method decreases exponentially. 

An essential aspect of any spectral method is the choice of 
expansion functions. Consider first the case of a bounded, cartesian 
domain. Fourier series are the most familiar expansion functions, but 
they are only appropriate for problems with periodic boundary 
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conditions. The appropriate collocation points on [ 0 , 2 tt ] are 

xj = 2Trj/N, j = 0,1,...,N-1. (7) 

In the general, non-periodic case, normalized to [-1,1], the appropriate 
class of functions is the Jacobi polynomials. The proper collocation 
points are generally -1, +1 and the extrema of the last polynomial 
retained (CHQZ, Ch. 2). The most commonly used Jacobi polynomials are 
the Chebyshev and Legendre ones. 

On an unbounded domain, the obvious choices of Laguerre or Hermite 
polynomials are rarely advisable. Not only are fast transforms unavail- 
able, but these expansion functions have relatively poor resolution 
properties (GO, Ch. 3). A better approach is to combine a mapping with 
a Fourier or Chebyshev series in the mapped variable. Boyd (1986) has 
shown that spectral accuracy can be achieved for u(x) on (-», ") 
with the mapping x = x A cot 5, and a full Fourier series in £, 
provided that u(x) exhibits at least algebraic decay at °°. Moreover, 
if u(x) has exponential decay, then a Fourier cosine series will 
suffice. The latter case is equivalent to a Chebyshev series in q with 
x = x* n//(l - n^). Spalart (1984) noted that the odd (or even) 
Chebyshev polynomials work well on [0,“), when combined with an expo- 
nential mapping, provided that u(x) decays faster than exponentially. 

The process of numerical differentiation is particularly simple 
when the expansion functions are trigonometric polynomials. Starting 
from u j , the values of u at x j , one computes 



7 


N-l 

a k = (1/N) ^ Uj exp(-ikx^ ) , k = (8) 

j=0 

and then uses 


N/2-1 



k=-N/2 


ik a^ exp ( ikx ) 


(9) 


to approximate du/dx at x j . The Fast Fourier Transform (FFT) can be 
used to evaluate both the sums given above. The total cost of computing 
the derivative in this manner is 5N log 2 N + N real operations. (All 
operation counts given in this review presume, for simplicity, that N 
is a power of 2 and that the complex FFT is used; however, FFTs which 
allow prime factors of 3 and 5 are just as efficient and are widely 
available and real to half-complex FFTs offer a 20% savings (Temperton 
(1983)).) The FFT can also be used to differentiate functions which are 
expanded in Chebyshev series, since expansions in these special Jacobi 
polynomials reduce to cosine series. Morever, in terms of the Chebyshev 
coefficients, derivatives are obtained by simple recursion relations 
(CHQZ, Ch. 2). For Chebyshev series the total operation count for 
differentiation is 5N log 2 N + 16 N. 

For the classical expansion functions the matrix which represents 
differentiation, i.e., d^u/dx^ = D^u, is known in closed form 
(Gottlieb, et al (1984)). Unlike the differentiation matrices for 
alternative, local discretizations, these matrices are full. Hence, the 
matrix-vector multiplication which produces the derivative at the collo- 
cation points costs 2N 4, operations. These operation counts suggest 
that for N > 16, transform methods are faster for differentiation 
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than matrix multiplies. On modern scalar and vector computers the 
transform methods become faster than the matrix-vector multiply methods 
for N between 16 and 32 (CHQZ, Ch. 2). 

An important issue in many applications of Chebyshev spectral 
methods is the manner in which the boundary conditions are enforced. 
Dirichlet boundary conditions are generally straightforward. Neumann 
boundary conditions may be enforced by altering the boundary values to 
ensure the desired normal derivative or by building the boundary 
condition into the differential operator (Streett, et al (1985)). For 
hyperbolic systems, characteristic boundary conditions are a virtual 
necessity (Gottlieb, et al (1981), Salas, et al (1985)). Canuto and 
Quarteroni (1986) discuss how to implement characteristic boundary 
conditions for implicit time discretizations. Chebyshev spectral 
methods have the advantages (over standard finite difference schemes) 
that they require the same number and type of boundary conditions as the 
analytical formulations of the problem, and that no special difference 
formulae are required at the boundary. 

The spectra of the discrete differentiation operators are an 
important characteristic of numerical methods. For Fourier approxima- 
tions to periodic problems, these are obvious: purely imaginary and 
growing as N/2 for D , negative real and growing as N^/4 for D . 
Indeed, for periodic problems such as u t + u x = 0, the Fourier 
eigenvalues are exactly equal to their analytic counterparts. This 
means that Fourier spectral methods propagate the numerical solution 
with zero phase error. This is illustrated in Figure 2 for the problem 
whose solution is u(x,t) = sin(u cos(x-t)). The lagging phase of the 
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finite difference solution is apparent, whereas the Fourier solution is 
indistinguishable from the true one. Of course, in realistic problems, 
variable coefficients or nonlinear terms will introduce non-zero (but 
still relatively small) phase errors. 

Figure 3 displays the eigenvalues of a Chebyshev approximation to 
d/dx on [-1,1] with homogeneous Dirichlet boundary conditions at 
x = +1. The eigenvalues are predominantly imaginary but do have 
negative real parts. The absolute value of the largest eigenvalue grows 
as N^. These eigenvalues may be surprising at first sight. After all, 
the periodic discrete problem has purely imaginary eigenvalues, whereas 
the non-periodic continuous problem has no discrete eignevalues. 
Nevertheless, Figure 3 does convey the nature of the eigenvalues of the 
discrete problem and these are crucial for both time differencing 
methods and iterative schemes. The eigenvalues of Chebyshev 
approximations to d^/dx^ with homogeneous Dirichlet boundary 
conditions at x = -1 and x = +1 are real and negative and the 
largest eigenvalue grows as (Gottlieb and Lustman (1983)). 

In practice, when one is solving an evolution problem such as u t = 
Lu, where the operator L contains all the spatial derivatives, one 
combines a spectral discretization of L with a standard finite 
difference technique for the time derivative. The Leap Frog, Adams- 
Bashforth, Crank-Nicolson, and Runge-Kutta schemes are the ones most 
commonly used (CHQZ, Ch. 4). The stability regions of these schemes 
depend upon the spatial operators. The stability properties of Fourier 
methods are qualitatively the same as those for second-order central 
difference spatial operators. However, the precise stability limit is 
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typically a factor of (l/ir) n smaller for Fourier approximations, 
where n is the order of the highest spatial derivative which appears 
in L. 

The stability properties of Chebyshev methods are more subtle. For 
example. Leap Frog is unconditionally unstable for advection problems, 
such as u fc + u x = 0, since the discrete eigenvalues of the spatial 
operator have negative real parts. On the other hand, second-order 
Adams-Bashforth and Runge-Kutta methods are strictly stable (and not 
weakly unstable like their Fourier counterparts) for the same reason. 

The typical time-step limitations on Chebyshev methods are 1/N^ 
for first derivative operators and 1/N^ for second derivative ones. 
These are far more stringent than the analogous restrictions for uniform 
grid finite difference approximations. They arise from the crowding of 
the collocation points near the boundaries (see Figure 1). Although 
this crowding necessitates small time-steps, it is required for the high 
spatial resolution of the method and is quite advantageous for problems 
with boundary layers. This is, however, a substantial disadvantage for 
problems with very little structure near the boundaries. It can be 
alleviated to some degree by mapping, but a mapping to a uniform grid is 
counter-productive because it destroys the spatial accuracy. 

This Chebyshev time-step limitation disappears when implicit time 
discretizations are employed. The principal difficulty is obtaining 
efficient solutions of the resulting implicit equations, since the 
matrices which represent the differentiation operators are full. In 
some special cases, fast direct solution methods are available. These 
typically require low-order polynomial coefficients and, in multi- 
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dimensional problems, at most one non-periodic direction (GO, Chs. 9 and 
10, Moser, et al (1983)). 

The use of implicit techniques in more general situations requires 
iterative methods. This has been one of the major developments of the 
current decade (CHQZ, Ch. 5). Let us denote a typical linear, implicit 
system arising from a spectral discretization as L g p u = f . The 
simplest iterative scheme — Richardson's method — is just 

u u + toff - L u) (10) 
v sp J 

where to is an acceleration parameter. The matrix L g p will be full 
and will have eigenvalues which grow rapidly as the number of grid 
points increases. The fullness of the matrix does not preclude 
iterative methods since transform techniques for differentiation permit 
the matrix-vector product L g p u to be computed in 0(N log2 N) 
operations rather than O(N^). The slow convergence which results from 
the large eigenvalues of L g p can ame li or-a t e< i by preconditioning. 
In this case the basic iterative scheme is 

u«-u + u>H^(f-L u) (11) 

where H is a preconditioning matrix. This will accelerate convergence 
if H is a good approximation to L g p, anc ^ be relatively 
inexpensive if H is readily inverted. The former condition is met by 
low-order finite difference (Orszag (1980)) and finite element (Deville 
and Mund (1985)) approximations to L g p. Although the latter condition 
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certainly holds for one dimensional problems, these particular 
preconditionings become increasingly expensive to invert as the 
dimensionality of the problem increases. 

The most attractive approach to very large problems is to combine a 
less accurate but more readily inverted preconditioning with multigrid 
techniques. Spectral multigrid methods take advantage of the fact that 
most iterative methods are highly effective in reducing the error 
components corresponding to the upper half of the eigenvalue spectrum, 
but are very inefficient on the remaining, low frequency components. 
Thus, in a multigrid method one combines iterations on the desired grid 
with (much cheaper) iterations on successively coarser grids. The 
details of this method are admittedly subtle, but they have been 
carefully described in a series of papers (Zang, et al (1982, 1984), 
Streett, et al (1985), Phillips, et al (1986)). Brandt, et al (1985) 
have demonstrated that many periodic problems can be successfully solved 
in this manner without the need for any preconditioning. 

Another recent innovation has been the development of spectral 
multidomain techniques. These allow spectral methods to be applied to 
geometries for which a single, global expansion is either impossible, or 
else inadvisable due to resolution requirements which vary widely over 
the domain. In a multidomain technique the full domain is partitioned 
into (not necessarily disjoint) subdomains. These may be patched 
together at interfaces or else they may overlap. The crucial part of 
the patched multidomain methods are the interface conditions. These may 
be expressed explicitly as continuity conditions (Orszag (1980), Kopriva 
(1986)), may arise from a variational principle (Patera (1984)), may 
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consist of integral constraints (Macaraeg and Streett (1986)), or may be 
enforced by a penalty method (Delves and Hall (1979)). The spectral 
element method of Patera is to date the most highly developed of 
these. Many techniques, such as isoparametric elements (Korczak and 
Patera (1986)), have been borrowed from conventional finite element 
methodology. Indeed, there are many similarities in this approach to 
the p-version of the finite element method (Babuska and Dorr (1981)). 
Figure 4 illustrates a spectral element grid as well as the computed 
solution for flow past a cylinder (Karniadakis , et al (1986)). In all 
cases convergence is achieved with a fixed number of subdomains while 
the number of grid points on each subdomain increases. The spectral 
overlapping subdomain methods were devised by Morchoisne, et al (1983), 
and are currently being investigated extensively in Europe. 


INVISCID FLOW 

Perhaps the simplest fluid dynamical problems are those which are 
steady, inviscid, incompressible and irrotational. In terms of the 
velocity potential $, these are described by the Laplace equation 

V 2 4> = 0 (12) 

with Neumann conditions on the boundaries. Spectral methods can be 
quite effective on such elliptic problems and also on the slightly more 
general class of problems described by 


V 2 <f> " X* = f 


(13) 
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with Dirichlet , Neumann or mixed boundary conditions. These more 
general methods could easily be applied to the idealized flow problem 
described above. 

Spectral methods have been developed for such Poisson/Helmholtz 
problems in a variety of geometries. Direct methods are straightforward 
when at most one of the directions requires non-periodic boundary 
conditions, and hence a Chebyshev polynomial representation. Constant 
coefficient equations become diagonal in the periodic directions. In a 
cartesian non-periodic direction the equation can be reduced to a quasi- 
l-^idi^gonal form (GO, Ch. 10) if the domain is finite and a penta— 
diagonal form if it is infinite and the cotangent mapping is used (Cain, 
et al (1984)). Otherwise, a matrix diagonalization technique can be 
employed (Murdock (1977), Haidvogel and Zang (1979)). Direct methods 
for problems with 2 or more non— periodic directions have been discussed 
by Haidvogel and Zang (1979), Haldenwang, et al (1984), and LeQuere and 
Roquefort (1985). Some extensions to three non— periodic directions are 
described by Haldenwang, et al (1984) and Tan (1985). Iterative methods 
allow efficient treatment of more general geometries, especially for 
exterior problems. See Canuto, et al (1985) and Deville and Mund (1985) 
for some standard techniques. Especially for very large problems of 
this type, spectral multigrid methods appear to be the most efficient 
(Zang, et al (1982, 1984)). 

Compressible potential flow is described by a similar, but 
nonlinear equation 


V • (p V <}>) = 0 


(14) 
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where the density p is a quadratic function of V<J> . For subsonic 
flow this problem is elliptic. Streett, et al (1985) have demonstrated 
the great efficiency that spectral multigrid methods achieve for this 
case. They have applied these techniques to the two-dimensional flow 
past a circular cylinder. Using a mere 2000 grid points they have 
obtained an estimate for the free stream Mach number at which the flow 
first becomes sonic. It agrees to 6 digits with the results of van Dyke 
and Guttman (1983) based on a Rayleigh-Jensen expansion. 

For transonic flow, the potential equation is of mixed type, with a 
supersonic pocket embedded in a subsonic flow. There will be a sonic 
line and usually a shock which terminates the supersonic region. The 
challenging numerical task is to obtain a converged solution to the 
discrete, nonlinear potential equation. Spectral multigrid methods have 
proven competitive with finite difference methods and achieve substan- 
tial economies in storage (Streett, et al (1985)). 

Still within the confines of inviscid flow, one can obtain the 
effects of vorticity by resorting to the Euler equations 

|f- + V • (pq) = 0 

9q. , 

_ + a . Vcj_= - f Vp. (15) 

|f + q • vs = o, 


where is the velocity, 
v s/ s n 

p = p 1 e . As is the 


p is the pressure, S is the entropy, and 
case for all numerical methods, the real 
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delicacy is the treatment of sonic lines and shock waves. The 
discontinuities arising from shocks are especially troublesome for 
spectral methods. The global nature of these approximations induces 
oscillations in the solution which are essentially of a Gibbs phenomenon 
type. The high frequency component of the solution decays very slowly. 
This part of the spectrum must be filtered to produce a presentable 
approximation. A detailed mathematical analysis of filtering techniques 
in Fourier spectral methods for linear, hyperbolic problems with 
discontinuous solutions has been presented by Majda, et al (1978). A 
post-processing procedure that involves matching the computed solution 
with simple discontinuities has been discussed by Abarbanel, et al 
(1985). 

The first applications of spectral methods to compressible flows 
focused on the treatment of shock waves in one-dimensional problems 
(Gottlieb, et al (1981), Zang and Hussaini (1981), Taylor, et al 
(1981)). As is the case with finite difference methods, spectral 
methods for problems involving shocks require some type of explicit or 
implicit numerical dissipation. In solutions to partial differential 
equations the explicit dissipation may take the form of a linear, 
spectral filter, or it may consist of an artificial viscosity term which 
is added to the Euler equations. This artificial viscosity may be 
nonlinear. Approximations based on Chebyshev polynomials may be stable 
without any explicit dissipation since the Chebyshev derivative operator 
contains implicit dissipation (Gottlieb, et al (1981)). 

Most investigations have confined themselves to problems whose 
solutions (even in two-dimensions, Sakell (1984)), were either piecewise 
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constant or else piecewise linear. No one has yet exhibited a spectral 
solution to a problem with both shock waves and complex flow structure 
in which spectral accuracy was attained (Hussaini, et al (1985a)). 

The difficulties that shock capturing spectral methods exhibit are 
not due to any intrinsic difficulty in resolving transonic and super- 
sonic flows. Kopriva, et al (1984) solved the Ringleb flow problem, 
which is a smooth two-dimensional transonic flow with a closed form 
solution, by a Chebyshev spectral method. They were able to exhibit the 
usual spectral accuracy on this class of problems. 

The shock-fitting approach popularized by Moretti (1968) for finite 
difference schemes was adapted to spectral discretizations by Salas, et 
al (1982). This technique avoids the Gibbs phenomenon by treating the 
shock as a boundary rather than as an interior region of the flow. It 
is applicable to flows which contain a few, geometrically simple shocks. 
Some problems for which high resolution results have been obtained by 
this method are the shock-vortex iteraction (Salas, et al (1982)), the 
shock-turbulence interaction (Zang, et al (1983), and the blunt body 
problem (Hussaini, et al (1985b)). 


BOUNDARY LAYER 

In many aerodynamic applications the boundary layer equations are 
an economical and useful model of viscous effects, especially when 
coupled interactively with an inviscid model for the outer flow (AGARD, 
1981). In similarity variables the two-dimensional boundary layer is 
described by 
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8 r 9 f . 3f , 2 3f 

3n 3n^ v 3n ^ 2?f 3£ “ 0 


3v 


3f 


IFf + f + 2? w - 0 


(.16) 


where f is the non-dimensional streamwise velocity, v is the normal 

velocity, n is the normal coordinate, and £ is the streamwise 

coordinate. The boundary conditions are f = v = 0 at n=0 and 
f + 1 as n + An inflow condition is required at some £ = £q* 

Chebyshev spectral approximations to the similar version of this 
system (3/35 = 0) are fairly straightforward to obtain by simple 
preconditioned iterative schemes (Streett, et al (1984)). This work 
demonstrated that a combination of domain truncation (typically at 
ri = 15) and grid stretching (to pack grid points near the solid 

boundary) is quite effective. A mere 20 collocation points will 

usually yield values for the wall shear and displacement thickness that 
have 3-digit accuracy, while 30 points produce 5-digit accuracy. 

The full non-similar equations are more challenging since there is 
a Chebyshev approximation in two directions. Streett, et al (1984) used 
an alternating direction type of preconditioning to obtain a solution. 
For non-similar flow roughly 20 polynomials in 5 (coupled with 25 in 
n) are required for 3-digit accuracy. They found that the Chebyshev 
approximation in 5 produced a substantial improvement over a simpler, 
mixed scheme which used finite differences in 5 together with 
Chebyshev collocation in p. The global nature of the streamwise 
approximation is especially useful for handling separated flow. In this 
case marching techniques in 5 are ineffective for finite difference 
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approximations. Figure 5 displays the streamlines and skin friction 
from a fully spectral solution of separated boundary layer flow. The 
arrow marks the region of the flow which is most sensitive to the 
numerical resolution. To obtain 4-digit accuracy in the skin friction 
here requires 40 collocation points in the normal direction and 26 in 
the streamwise direction. The corresponding requirements for a standard 
second-order finite difference method are 240 and 200 points 
respectively. Moreover, the spectral solution requires only 10% of the 
CPU time taken by the finite difference method. 


NAVIER-STOKES FLOW 

Much of the current enthusiasm for spectral methods is attributable 
to their success on simple, yet computationally intensive problems in 
viscous, time-dependent, incompressible flow. The pioneering 
simulations of three-dimensional homogeneous, isotropic turbulence by 
Orszag and Patterson (1972) were particularly influential. Subsequent 
calculations of three-dimensional transition and turbulence in simple 
wall-bounded flows have also been persuasive. Algorithms for these 
problems are substantially more difficult and time-consuming than those 
for homogeneous flows. The presence of non-periodic boundary conditions 
makes purely Fourier methods inappropriate and detailed simulations of 
transition problems typically require an order of magnitude more time- 
steps than do turbulence problems. The simplest class of such problems 
consists of flows which are assumed to be periodic in two directions, 
e.g., Poiseuille flow and Taylor-Couette flow for cylinders of infinite 
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length. In these cases, ones needs a Chebyshev discretization in only 
one direction. Several types of multidomain spectral techniques are 
currently being explored to extend further the class of viscous problems 
that are amenable to spectral methods. 

In many applications the preferred version of the Navier-Stokes 
equations is 


8q 

Jt 


U x £ = - VP + vv' 


<L 


(17) 


= 0 . 


The velocity is denoted by _c[_, the pressure by p, the vorticity by 
(o = V x u, the pressure head by P = p + (1/2(%| 2 , and the kinematic 
viscosity by v. This so-called "rotation" form is favored because, as 
noted by Orszag (1972), the use of the rotation form guarantees that 
Fourier collocation methods conserve kinetic energy. One can easily 
show that momentum is conserved as well. The conservation of kinetic 
energy is especially important for numerical reasons. In practice, it 
means that if the time differencing scheme is operated at time-steps 
below its stability limit, then nonlinear instabilities will not occur. 


Homogeneous Turbulence 

Homogeneous, isotropic turbulence is perhaps the one fluid 
dynamical problem for which strictly periodic boundary conditions in all 
spatial directions are justifiable. Hence, Fourier spectral methods are 
ideally suited for this class of problems. Moreover, since the 
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nonlinearities of the Navier Stokes equations are at worst quadratic, 
Fourier Galerkin methods are the most natural and efficient spectral 
technique for this problem (Orszag and Patterson (1972). Rogallo (1981) 
developed a linear coordinate transformation that permits simulation of 
flows with constant strain, shear and rotation within the confines of 
periodic boundary conditions. Rogallo (1981) and Basdevant (1983) have 
discussed techniques for minimizing the storage, CPU time and I/O costs 
of such algorithms. 

3 

The original simulations of Orszag and Patterson were on 32 
grids. By the early 1980's, 64 3 simulations were fairly routine. 

Rogallo (1981), Kerr (1985) and Lee and Reynolds (1985) have performed 

3 

numerous 128 simulations. By fully exploiting the special symmetries 
of the Taylor-Green vortex, Brachet, et al (1983) achieved a simulation 
of this flow at a Reynolds number of 3000 with an effective resolution 
of 256 3 . 

Fourier collocation approximations to this problem are also 
possible. For these approximations use of the rotation form of the 
Navier-Stokes equations is crucial. (Galerkin approximations to the 
inviscid form of these equations will automatically conserve momentum 
and kinetic energy in the absence of time differencing errors). 

The review by Rogallo and Moin (1984) discusses many applications 
of these techniques to problems in homogeneous turbulence. Here we need 
mention only the most recent applications. A primary goal of most of 
the simulations of isotropic turbulence has been to establish 
numerically the existence of an inertial range. The inertial range has, 
of course, been established experimentally, but only for Reynolds 
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numbers exceeding 10,000. Even though the high resolution calculations 
of Brachet, et al (1983) were performed at a Reynolds number of 3000, 
which is uncomfortably low by experimental standards, they did achieve 
the first plausible inertial range in a numerical simulation of 
turbulence. Bardina, et al (1985) and Dang and Roy (1985) have 
simulated the evolution of turbulence intensity in rotating flow. Wu, 
et al (1985) have performed calculations of compressed turbulence. Lee 
and Reynolds (1985) have analyzed the structure of turbulence in 
axisymmetrically contracting and expanding flow. Moin, et al (1985) 
have used numerical simulations to extract the large-scale vortical 
structures of some turbulent shear flows. Kerr (1985) has examined 
high-order correlations and small-scale structure in isotropic 
turbulence involving passive scalars. 

A few applications, all using the collocation technique, have been 
made to compressible, homogeneous turbulence. Feiereisen, et al (1981) 
simulated subsonic turbulent flows with uniform shear. They used a 
collocation method, in part because a Galerkin method is much more 
cumbersome and costly for problems with more than quadratic 
nonlinearities. Compressible, two-dimensional turbulence has been 
investigated by Leorat, et al (1985) and by Delorme (1984), the former 
with a fairly standard scheme and the latter with an implicit time- 
differencing method based on the ideas of Lerat, et al (1982). 

Linear Stability 

Most investigations of stability and transition in wall-bounded 
flows rely, at least in part, upon the results of linear stability 
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theory. The Orr-Sommerfeld equation has been the basis for many 
investigations of the stability of incompressible parallel flows (Drazin 
and Reid (1981)). This eigenvalue problem is described by a fourth- 
order ordinary differential equation. The Chebyshev approximation 
developed by Orszag (1971) for the temporal stability problem has been 
adopted and extended by many investigators. (A separate development of 
Chebyshev methods for ordinary differential eigenvalue problems has been 
conducted by Ortiz. See Chaves and Ortiz (1968) and, more recently, 
Ortiz and Samara (1983).) Leonard and Wray (1982) developed a Galerkin 
method for pipe flow which uses special Jacobi polynomials. Spalart 
(1984) demonstrated that for exterior flows (such as the parallel 
boundary layer) the use of only half the usual Chebyshev basis was 
advisable. Boyd (1985) has developed methods in the complex plane which 
are useful for flows in which the critical layer is well-separated from 
the wall. Von Kerczek (1982) has used Chebyshev polynomials for 
assessing the stability of oscillatory plane Poiseuille flow. Mac 
Giolla Mhuiris (1986) has used the Galerkin technique to examine the 
linear stability of some axisymmetric flows which are relevant to the 
vortex breakdown problem. 

The spatial stability versions of these problems are more difficult 
because the eigenvalue enters nonlinearly. Chebyshev methods for time- 
independent but spatially growing perturbations of Poiseuille flow are 
discussed by Bramley and Dennis (1982). Bridges and Morris (1984a, 
1984b) solved by a spectral method the more difficult, general spatial 
stability problem of self-similar boundary layers. 

These methods have been extended, in the manner of Floquet theory, 
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to include weakly nonlinear effects. In addition to a Chebyshev 
discretization in the direction normal to the wall, one includes several 
Fourier harmonics in the streamwise direction. Orszag and Patera (1983) 
and Herbert (1983a) have used this approach to determine the neutral 
stability surface of finite amplitude two-dimensional Tollmien- 
Schlichting waves in channel flow. In turn, the linear stability of 
these neutral finite amplitude waves can be examined. Thus, the linear 
stability of some special, temporally and spatially varying flows can be 
examined. Orszag and Patera (1983) have used this technique to study 
the interaction of two-dimensional and three-dimensional Tollmien- 
Schlichting waves in channel flow. Herbert (1983a, 1983b, 1984) has 
performed a detailed study of channel and boundary layer flows. He has 
unravelled the details of fundamental and subharmonic instabilities in 
parallel flows. 

Transition 

Transition to turbulence is highly nonlinear and a full simulation 
of the Navier-Stokes equations is required for its investigation. The 
primary difficulty of algorithms for incompressible flows is the 
simultaneous enforcement of the incompressiblity constraint and the no- 
slip boundary condition. This constraint is most easily but least 
rigorously satisfied in splitting methods, of which the Orszag-Kells 
(1980) algorithm is the prototype. The splitting errors of this method 
are 0(1) near the boundary for the normal pressure gradient and 
diffusion terms (Deville 1985). They appear to cause no serious errors 
in the channel flow problem. However, Marcus (1984a) decisively 
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demonstrated that the boundary errors produce serious inaccuracies in 
Taylor-Couette flow — as both the spatial and temporal discretizations 
are refined, the algorithm appears to converge to answers that disagree 
with experiments in the third digit. Marcus (1984a) and Kleiser and 
Schumann (1984) devised an influence matrix technique which completely 
eliminates the splitting errors at a modest extra cost. Marcus found 
that the results of this algorithm agreed with the experimental results 
to the full 4 digits that were available. He ascribed the sensitivity 
of the rotating cylinder problem to the fact that its dynamics are 
driven by the motion of the boundary rather than by a mean pressure 
gradient. 

A procedure for reducing, although not entirely eliminating the 
splitting errors at the boundary, was devised by Fortin, et al (1971) 
for finite element methods (re-discovered later by Kim and Moin (1985)), 
and applied to spectral algorithms by Zang and Hussaini (1986). It 
consists of modifying the boundary conditions for the intermediate steps 
of the algorithm so that both the no-slip and divergence-free conditions 
are satisfied at the end of the full time-step to a higher order in the 
size of the time-step. 

The big advantage of these splitting techniques is that they 
require the solution of only Poisson equations (for the pressure) or 
Helmholtz equations (from a Crank-Nicolson discretization of the viscous 
term). These positive definite, scalar equations are much easier to 
solve numerically that the indefinite, coupled equations that arise in 
unsplit methods. The Orszag-Kells , Marcus and Kleiser-Schumann 
algorithms resort to direct solution methods of the type discussed in 
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section 3. The Zang— Hussaini algorithm employs iterative techniques so 
that it is applicable to a wider class of problems. The most sophis- 
ticated and powerful of the iterative techniques is the spectral multi- 
grid method. It makes the cost of a single time— step of order 

3 

N log 2 N, even for problems with variable geometric terms and transport 
coefficients. In contrast, a parallel flow problem, even with uniform 
transport coefficients, requires order operations per step by 

direct methods. 

One way to avoid the splitting errors is to integrate the 
incompressible Navier-Stokes equations in a single step that couples the 
divergence-free constraint with the momentum equations. The numerical 
difficulty of this approach is that one must invert a larger set of 

equations (it involves the pressure as well as the three velocity 

components), which is indefinite. In a few special cases direct 

techniques are viable (Moin and Kim (1980)). The preconditioned 
iterative scheme of Malik, et al (1985) has been applied to channel flow 
(Zang and Hussaini (1985a)) and to the heated boundary layer (Zang and 
Hussaini (1985b)), a problem which involves variable transport 
coefficients, and also in a verification of weakly , nonlinear stability 
theory for stagnation point flow (Hall and Malik (1986)). 

Many of the numerical problems caused by the incompressibility 
constraint can be avoided by an expansion in functions which are 
divergence-free (Ladyzhenskaya 1969, Temam 1977). Leonard and Wray 
(1982) first applied this idea to spectral methods. They devised a set 
of basis functions for pipe flow which are both divergence— free and 
satisfy no— slip boundary conditions. Similar basis functions have been 



27 


developed for straight and curved channels (Moser, et al (1983)) and for 
the parallel boundary layer (Spalart (1984)). This class of methods can 
be quite economical of storage since only two variables per grid point 
are required to specify the flow field. (However, in actual 

implementations it may be more efficient in terms of CPU time to store 
several additional quantities per grid point.) The efficiency of these 
methods is dependent upon the bandwidth of the matrices which arise from 
the implicit treatment of the viscous terms. In the examples cited 
above, the bandwidth is quite small, roughly of order 10. This 
requirement has dictated the use of special Jacobi polynomials rather 
than Chebyshev ones in pipe and boundary layer flow. As a consequence, 
transform methods are not applicable in the non-periodic direction. 
Hence, the cost of evaluating the nonlinear terms increases as 
rather than as l°g 2 N. Moreover, in even slightly more general 

cases, the matrices can be completely full. 

Orszag and Patera (1983) performed a parametric study of the 
secondary instability in channels and pipes, demonstrating that 
subcritical instabilities exist at Reynolds numbers as low as 1100. 
Kleiser and Schumann (1984) replicated many of the features of the 
Nishioka, et al (1980) experiments on channel flow transition. Both 

groups also obtained good quantitative agreement with the predictions of 
weakly nonlinear theory. The subharmonic instabilities that were 
predicted by Herbert's (1983b), (1984) weakly nonlinear analysis (and 

are also in evidence in boundary layer experiments, Saric, et al 1984), 
were reproduced by Spalart (1985) and Laurien (1986) for the boundary 
layer and by Zang and Hussaini (1985a), and by Singer, et al (1986) for 
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channel flow. The existence of a similar nonlinear instability of 
center modes in channel flow was uncovered by Zang and Hussaini 
(1985a). A detailed comparision of nonlinear effects on the laminar 
flow control techniques of pressure gradient, suction and heating in 
boundary layer flow was made by Zang and Hussaini (1985b). Krist and 
Zang (1986) have performed a detailed study of the resolution require- 
ments for simulation of the later stages of transition to turbulence in 
channel flow. The spanwise direction places the greatest demands on the 
resolution because of the very sharp spanwise gradients which occur near 
the tip of the characteristic hairpin vortex. Figure 6, which is 
extracted from that work, illustrates the structure. 

Marcus (1984a, 1984b) has performed a careful numerical study of 
non-axisymmetric instabilities in classical Taylor-Couette flow. He has 
produced 4-digit agreement with the wave speeds measured by King, et al 
(1984) for both the one wavy-vortex and the two wavy vortex states. 
Marcus and Tuckerman (1986a, 1986b) have simulated axisymmetric 
spherical Couette flow. Unlike previous workers they did not assume 
equatorial symmetry. This was a crucial factor in their success in 
reproducing the transitions between 0, 1, and 2 vortex states observed 
by Wimmer (1976). 

Inhomogeneous Turbulence 

In several cases these algorithms have been used to simulate 
turbulence in wall-bounded flows. Orszag and Patera (1983) performed a 
64 simulation of turbulent channel flow which reproduced the turbulent 
velocity profile, including the law of the wall behavior. Moser, et al 
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A 

(1984) computed turbulent flow in a curved channel on a 128^ x 64 

grid. They reproduced some of the data on low-order turbulence 
statistics and exhibited some of the effects of curvature. Spalart and 
Leonard (1985) have done some analyses of pressure gradient effects in 
turbulent boundary layers. 

More Realistic Geometries 

As noted above, there is a substantial increase in cost when there 
is more than one inhomogeneous direction in the problem. The Kleiser- 

Schumann influence matrix technique has been extended to two non- 

periodic directions by LeQuere and Roquefort (1985), who used it to 
study thermal convection in a square cavity. Streett and Hussaini 
(1986) similarly extended the split algorithm of Zang and Hussaini 

(1986), and used it to study the effect of finite length cylinders in 
Taylor-Couette flow. Ku and Taylor (1985) have developed an algorithm 
for three non-periodic directions. This method presently treats only 
the pressure term implicitly. Thus, there can be a severe time-step 
limitation arising from the viscous terms. Morchoisne (1984) has 

developed a number of methods for problems with more than one non- 
periodic direction. In general iterative techniques are used for 
solving the resulting implicit equations. There has not yet been any 
systematic comparison of these methods. Leonard (1984) has derived a 
set of divergence— free basis functions for 2 nonperiodic directions, but 
an efficient solution technique for the implicit equations has not yet 
been devised. 

Several of the multidomain spectral methods have been applied to 
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viscous problems. Morchoisne (1984) has performed some sample calcula- 
tions of channel flow. The spectral element has been used to calculate 
heat transfer in a two-dimensional, grooved channel (Ghaddar, et al 
(1984)) and to investigate stability and resonance phenomena in imbedded 
cavities in channel flows (Ghaddar, et al (1986a, 1986b)). Other 
applications include two-dimensional flow past a cylinder and flow past 
three-dimensional roughness elements (Karniadakis, et al (1986)). 

Spectral/Finite Difference and Quasi-Spectral Methods 

Heretofore, this review has been confined to numerical fluid 
dynamical work which employed spectral discretizations in all coordinate 
directions. There have, of course, been numerous computations which 
used mixed spectral/finite difference methods, i.e., algorithms with 
spectral discretizations in some directions and finite differences in 
the others. The parallel boundary layer transition calculations of Wray 
and Hussaini (1984) fall into this category. They used a Fourier 
spectral method in two periodic directions and second-order finite 
differences in the normal direction. They demonstrated that, despite 
the neglect of non-parallel effects, these simulations could reproduce 
features observed experimentally by Kovasznay, et al (1962), up to the 
so-called "two-spike stage" of transition. A slightly different 
spectral/finite difference method was used by Moin and Kim (1982) in 
their large-eddy simulations of turbulent channel flow and by Biringen 
(1985) in a study of active control in channel flows. More recently, 
Eidson, et al (1986) have used a similar algorithm in a high resolution 
direct simulation of a turbulent Rayleigh-Benard flow. 

Another alternative to true spectral methods are what might be 
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termed quasi-spectral methods. Such algorithms employ Fourier expan- 
sions in all directions but infinite-order accuracy is not attained due 
to non-periodic physical boundary conditions in at least one 
direction. The simulations by Riley and Metcalfe (1980) of a time- 
developing mixing layer fall into this category. In this idealized flow 
the mean velocity is solely a function of the transverse coordinate 
y. Although the flow extends to y = ±», Riley and Metcalfe computed on 
a finite domain in y and used sine or cosine expansion to enforce 
free-slip boundary conditions in y. Quasi-spectral methods have also 
been used by Curry, et al (1984) to study Benard convection. 

True spectral methods have been developed for the time developing 
mixing layer. Cain, et al (1984) use a cotangent transformation in y 
combined with a Fourier method. Metcalfe, et al (1986) use hyperbolic 
tangent or algebraic transformations combined with a Chebyshev method. 

Riley and Metcalfe (1980) find that large amplitude two-dimensional 
disturbances have a pronounced effect upon the evolution of a turbulent 
mixing layer. Metcalfe, et al (1986) find that the mixing layer 
exhibits three-dimensional secondary instabilities similar to those 
which occur in wall-bounded flows. They appear to account for the 
mushroom-shaped features which are observed experimentally. Cain, et al 
(1981) have performed large-eddy simulations of this problem. 


REACTING FLOWS 

An emerging application field for spectral methods is reacting 
flows. These flows are especially challenging because they contain 
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sharp gradients in both space and time and because most real flows 
involve dozens or even hundreds of species. Flame fronts and shock 
waves are an additional complication. Some of the important features 
are mixing rates, ignition, and flame holding. 

There are a number of simplifying assumptions which lead to more 
tractable, but less realistic, models of reacting flows. The most 
drastic of these is that the reactions proceed without heat release and 
that the Mach number is so low that the flow may be trated as incompres- 
sible. Riley, et al (1986) have performed some three-dimensional 

simulations of a 2 species, time-developing mixing layer. They used a 
quasi-spectral method and obtained good agreement with both similarity 
theory and experimental data. 

McMurtry, et al (1986) employed a low Mach number approximation 
which includes some mild heat release effects but neglects the acoustic 
modes. They performed some two-dimensional calculations which indicate 
that the first-order effect of heat release is to reduce the rate of 
mixing. 

Drummond, et al (1986) applied a Chebyshev spectral method to a 
supersonic quasi-one-dimensional diverging nozzle flow with a simple but 
quite stiff 2 species hydrogen-air reaction. The spectral method 
proved to be quite economical compared with a benchmark finite 
difference result. The Chebyshev grid point distribution was quite 
well-adapted to the sharp gradients at the nozzle inflow, but less well- 
suited to the fairly uniform outflow region. 



33 


PERSPECTIVE 

A decade ago spectral methods appeared to be well-suited only to 
problems governed by ordinary differential equations or by partial 
differential equations with periodic boundary conditions. And, of 
course, the solution itself needed to be smooth. Some of the obstacles 
to wider application of spectral methods were: 1) sensitivity to 
boundary conditions; 2) treatment of discontinuous solutions; 3) 
resolution and time-step limitations imposed by the standard spectral 
grids; and 4) drastic geometric constraints. 

Substantial progress has been made on the implementation of Neumann 
boundary conditions, on characteristic boundary conditions for 
hyperbolic systems, and on the use of pressure and intermediate boundary 
conditions in incompressible flow. There have been some theoretical 
advances on filtering techniques for discontinuous solutions to linear 
problems. Moreover, the development of shock-fitting techniques has 
opened a new field of applications to compressible flows with shock 
waves. Some efficient direct solution techniques have been devised 
which enable severe viscous time-step limitations to be overcome in 
certain special geometries. The development of preconditioned iterative 
methods and, in particular, spectral multigrid techniques have radically 
expanded the class of problems which can be handled efficiently by 
spectral methods. Moreover, they lend much greater flexibility 
(combined with mapping techniques) to the grid point distribution. 
Finally, various multidomain techniques have expanded the range of 
spectral methods to many problems of real, practical interest. 
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CAPTIONS 

1. Comparison of finite difference (left) and Chebyshev spectral (right) 
differentiation. The solid curves represent the exact function and the 
dashed curves their numerical approximations. The solid lines are the 
exact tangents at x = 0 and the dashed lines the approximate tangents. 
The error in slope is noted as is the number of intervals N. 


2. Finite difference (circles) and Fourier spectral (diamonds) approximations 
after one period to a simple wave equation whose exact solution is 
represented by the curve. 

3. Eigenvalues of the Chebyshev first derivative operator for N = 16. 

4. A. spectral element grid (top) and the corresponding numerical solution for 
flow past a circular cylinder (courtesy of G. E. Karniadakis and A. T. 
Patera) . 

5. Streamlines (top) and skin friction (bottom) from a Chebyshev spectral 
solution of the boundary layer equations (courtesy of C. Streett). 

6. Streamwise (left) and spanwise (right) vorticity at four streamwise 
locations for a hairpin vortex in low Reynolds number channel flow 
transition. Only the lower half of the channel is shown. 
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